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The problem of energy conservation in the lattice Boltzmann method is solved. A novel model 
with energy conservation is derived from Boltzmann's kinetic theory. It is demonstrated that the 
full thermo-hydrodynamics pertinent to the Boltzmann equation is recovered in the domain where 
variations around the reference temperature are small. Simulation of a Poiseuille micro-flow is per- 
formed in a quantitative agreement with exact results for low and moderate Knudsen numbers. The 
new model extends in a natural way the standard lattice Boltzmann method to a thermodynamically 
consistent simulation tool for nearly-incompressible flows. 
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The overwhelming majority of fluid flows of physical 
and engineering interest are slow. That is, characteristic 
flow speed u is small compared to the speed of sound c s . 
This is quantified by the Mach number, Ma ~ u/c B . In 
typical situations, Ma varies from 1CT 3 - 1CT 2 in hydro- 
dynamic flows (turbines, reactors etc) to 10~ 4 in flows 
at a micrometer scale. In this paper we address a wide 
class of flows at Ma <C 1. Then the simplest charac- 
terization of the degree of molecularity is the Knudsen 
number Kn ~ X/H, the ratio of the mean free path of 
molecules A and the characteristic scale H of variation of 
hydrodynamic fields (density, momentum, and energy). 
When Kn < 10 -3 , one considers the hydrodynamic limit 
where molecularity reduces to specific for each molecu- 
lar model set of transport coefficients (viscosity, thermal 
conductivity etc). If, in addition, the Mach number is 
also small, one enjoys the incompressible hydrodynamics 
with the ordering Kn <C Ma <C 1, and the flow can be 
characterized solely by the ratio Re ~ Ma/Kn (one of 
the definitions of the Reynolds number). 

Contemporary computational fluid dynamics becomes 
increasingly more interested in the domain where Mach 
number remains small but Knudsen number increases, 
thus, the incompressibility becomes gradually lost. Be- 
cause of its relevance to the engineering of micro electro- 
mechanical systems (MEMS), the branch of computa- 
tional fluid dynamics focused on micro scale phenomena 
is often called "micro-fluidics" . Typical flows in micro- 
devices are highly subsonic (with characteristic flow ve- 
locities about 0.2 m/s, corresponding to Ma ~ 10 -4 ), 
while Knudsen number varies from Kn ~ 10~ 2 (so-called 
slip-flow regime) to Kn ~ 1 (moderately rarefied gas 
flows) 0. There is much need for computational models 
in the domain of slow flows where the effects of molecu- 
larity become increasingly more pronounced. 

In recent years, the lattice Boltzmann method has 
drawn considerable attention as a simulation method for 
flows at low Mach numbers. Especially popular are the 
so-called isothermal lattice Boltzmann models (ILBM) 
without energy conservation 0]. Recently, there was in- 
creasing interest in applying these models also for micro- 



flow simulations 0, H 0, H, IE 00 The hydrodynamic 
(locally conserved) fields in the ILBM are the density 
p and the momentum density j, whereas the conserva- 
tion of the energy is not addressed. Construction of the 
ILBM may vary among the authors, but all these models 
have one point in common: The lack of energy conser- 
vation inevitably leads to a bulk viscosity. Indeed, the 
non-equilibrium part of the stress tensor in ILBM reads: 
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This tensor is not traceless, P£° q ~ 2Kn<9 Q (j Q / p), which 
immediately leads to the bulk viscosity terms in the equa- 
tion for the momentum density. We remind that the 
physical bulk viscosity in hydrodynamic models is related 
to a redistribution of the energy among the translations! 
and internal degrees of freedom of molecules rather than 
to any non-conservation of the energy. Thus, from the 
physical standpoint, the bulk viscosity of the ILB mod- 
els is spurious. Certainly, the presence of the bulk vis- 
cosity, spurious or not, by no means precludes the limit 
of incompressible hydrodynamics because, loosely speak- 
ing, the divergence of the velocity field u = j/p vanishes 
in that limit |9(. Thus, ILBM is a valid model for the 
incompressible hydrodynamics. However, the spurious 
bulk viscosity of ILBM becomes a severe drawback when 
such models are applied to weakly-compressible or micro- 
flow simulations. 

The best way to illustrate this problem is to consider 
a representative example. Plane Poiseuille flow is one of 
the most studied benchmarks on gas dynamics. The gas 
moves between two parallel plates driven by a fixed pres- 
sure difference between the inlet and outlet. It is known 
from the classical kinetic theory 01 that the flow rate 
Q has the following asymptotic at low and high Knudsen 
numbers: 



Qo = (6Kn) _1 + s + (2s 2 - I) Kn, Kn < I 
Qoo ~ (l/V^F)ln(Kn)+0(l), Kn->oo, 
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with s = 1.015. These two asymptotic limits ensure that 
the flow rate has a minimum at some finite Kn (the Knud- 
sen minimum [13 )• While a qualitative agreement of 
the ILBM simulations with the continuous-velocity ki- 
netic theory was found at all Knudsen numbers 0, j 
the quantitative agreement is poor beyond the slip-flow 
regime at Kn > 10~ 2 . It was found that the ILBM sys- 
tematically over-predicts the flow rate at small Knudsen 
numbers. This is the effect of the bulk viscosity which 
can be qualitatively explained as follows: At low Knudsen 
numbers the behavior is still dominated by particle's col- 
lisions in the bulk, therefore, the steady state is reached 
upon a balance between the frictional force ~ Knd/sP^p 
and the forcing due to the constant pressure gap between 
the inlet and the outlet. Therefore, if there is additional 
contribution of the bulk viscosity (more friction), this 
balance at the same Kn shifts to a higher velocity at the 
steady state, and will result in the over-prediction of the 
flow rate. On the other hand, the lack of the energy con- 
servation also contributes to the over-prediction at high 
Knudsen numbers by a different mechanism: Since colli- 
sions in the bulk become rear, it becomes important that 
the energy be correctly redistributed in these rear events. 

In this paper we introduce new lattice Boltzmann mod- 
els with the energy conservation. These models are de- 
rived from the continuous kinetic theory, are free from the 
drawbacks of the isothermal lattice Boltzmann models, 
and at the same time they retain in full the outstanding 
computational efficiency of the latter. 

The structure of the paper is as follows: First, we shall 
derive the discrete-velocity model from the continuous- 
velocity kinetic theory. Second, we shall redo the compu- 
tation of the micro-Poiseuille flow with the new model, 
and demonstrate a significant (roughly, order of magni- 
tude in Kn) improvement in the accuracy with respect to 
isothermal simulation. Third, we shall explain in which 
points the present derivation differs from the earlier stud- 
ies. A brief discussion concludes the paper. 

Starting point of our derivation is the grand canonic 
potential of the Boltzmann kinetic theory, 



H = J FlnFdv + p J Fdv + ( a J Fv a dv + j J Fv 2 dv, 
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where F(x,v) is the one-particle distribution function, 
and p, ( a , and 7 are Lagrange multipliers correspond- 
ing to density, momentum and energy, respectively. The 
D + 2-parametric family of functionals J5J, where D is 
the dimension of the velocity space, describes the equi- 
librium states as its minima, 8H — 0, and it also defines 
the locally conserved fields (density p, momentum j, and 
energy e), 
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In order to derive the discrete velocity kinetic theory. 



the functional 101 is evaluated with the help of the D- 
dimensional Gauss-Hermite quadrature with the Gaus- 
sian weight W = (2n9 Q ) 2 exp (— (mv 2 )/(26> )), where 
$0 = {kBTo/m) is the reduced uniform reference tem- 
perature. We remind that the quadrature evaluation 
of an integral replaces it by a sum, J VK(v)G(v)c?v w 
Y^,7=i WiG(yi), where Vj, i= 1, . . . , na are the nodes (or 
abscissas) of the quadrature, and W% are corresponding 
weights. In the case under consideration, the nodes of 
the quadrature (discrete velocities) are situated at the 
zeroes of Hermite polynomials. 

We now proceed with the Gauss-Hermite quadrature 
evaluation of the velocity integral J5J. For concrete- 
ness, we shall consider the third-order Hermite polyno- 
mial. Then rid = 3 s , and the discrete velocities and 
weights is constructed as follows: For D = 1, the three 
roots and corresponding weights are (—^/39 , 0, y/30 ), 
(1/6, 2/3, 1/6); for D > 1, the roots are all possible ten- 
sor products of the roots in D = 1, and the weights 
are corresponding products of one-dimensional weights. 
We shall consider D = 3 below, that is rid = 27 (same 
considerations apply to any quadrature, in particular, to 
the popular 9-velocity model for D = 2). As is well 
known, the third-order quadrature has the unique feature 
that its nodes form a face-centered square lattice which 
is the crucial feature to the further lattice Boltzmann 
discretization in space and time. Introducing the pop- 
ulations, fi = W^(27r6» ) 3 / 2 exp(w 2 /(26l ))^(a;,v 4 ), and 
using the reduced discrete velocities, Cj = Vj/-\/30o, we 
write the quadrature for J3J) 
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Differentiation of (@J with respect to Lagrange multipliers 
defines the locally conserved fields in the discrete case (cf. 
©), 
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XX 1 ' Cm ' = { p > 3p + P~ 1 j 2 }- (5) 
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The equilibria /° q are now found as minima of H 
From the extremum condition. 5H = 0, it follows 



= W i ex P{-^ - (aC ia - 7C 2 }. 
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In order to express the Lagrange multipliers in JSJ) in 
terms of hydrodynamic fields J5J), we substitute © into 
(J5J and derive the functions p,{p,j,p), Ca(piJ)P) an d 
l{Pi3iP) by perturbation for small momentum, owing 
for the fact that Ca(p,0,p) = 0, and that the zero- 
momentum functions p(p,0,p) and j(p,0,p) can be 
found in a closed form. Computation is quite straightfor- 
ward, and we write here the final result to second order 
in the momentum: 
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The pre-factor in this formula has the following limit 
when (p/p) — * (1/3): 
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The implication of this limit will be important below 
when we will be discussing the relation of the present 
model to the isothermal lattice Boltzmann model. 

We now proceed with the evaluation of the stress ten- 
sor P^piPiJiP) an d °f the energy flux q^(p,j,p) at equi- 
librium. The important observation to be made here is 
that if the pressure to density ratio satisfies the condi- 
tion, then P^p and satisfy the corresponding rela- 
tions pertinent to the continuous-velocity Maxwell dis- 
tribution. In dimensional units, the condition just men- 
tioned reads p = (k&Top)/m, that is, it corresponds to 
the ideal gas equation of state at the reference tempera- 
ture of the Gaussian weight. Moreover, if we allow small 
variations of the pressure around the point p/p = 1/3, 



namely, \{p/p) - (1/3) j 

a/3 



Ma , the Maxwell's form of the 



functions P®% and q® q persists, and we have 
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Note that the variation of the pressure to density ratio 
of the order Ma 2 is pertinent to flow phenomena where 
the dominant effect on the temperature is the viscous 
dissipation and associated heat conduction. Condition 
I (p/p) ~ (1/3)1 ~ Ma 2 is a conservative estimate of the 
domain of validity of the present model. 

With the equilibrium J7J, we write up the simplest 
kinetic equation (the Bhatnagar-Gross-Krook model) , 

dth + dadcfi = ~-(fi ~ fP(p,j, P )), (11) 
T 

where r > is the relaxation time. In order to find out 
the hydrodynamic limit of the model, we perform the 
Chapman-Enskog analysis at low Mach numbers. In so 
doing, we neglect all terms in j a of the order three and 
higher, and we end up with the following non-equilibrium 
expressions for the stress and the heat flux: 




1/Kn 



FIG. 1: Flow rate in the pressure driven Poiseuille flow as 
a function of inverse Knudsen number. Comparison of the 
present energy-conserving model with the isothermal lattice 
Boltzmann model D2Q9 and the continuous linearized 
Boltzmann-BGK model 
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The most important achievement is that the non- 
equilibrium (Newtonian) stress (|12f) is traceless, that is, 
by preserving the energy conservation we eliminated the 
spurious bulk viscosity. The heat flux l|13J) obeys the 
Fourier law. 



We have implemented the lattice Boltzmann space- 
time discretization of the kinetic equation i|ll|) , and redo 
the micro-Poiseuille flow simulation mentioned in the in- 
troduction. Results are presented in Fig. ??, where the 
present model is compared to the exact solution of the 
continuous linearized BGK model [l4|, and the 2DQ9 
isothermal lattice Boltzmann model with the spurious 
bulk viscosity 13]. It is clearly visible that the effect 
of the spurious bulk viscosity of the isothermal model is 
completely eliminated, and that the quantitative agree- 
ment with the continuous BGK model extends up to 
Kn = 0.1. 
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Finally, let us place our derivation with respect to pre- 
viously reported lattice Boltzmann models on the same 
lattice. If we substitute p — (l/3)p into the equilibrium 
function Q, and use the limit ©, then / 4 eq (p, j, p/3) 
recovers the second-order polynomial equilibrium of the 
isothermal lattice Boltzmann method on the same lat- 
tice 0|, and instead of the traceless stress tensor 
()12J) we recover Q with the bulk viscosity component. It 
needs to be stressed that the second-order polynomial in 



j (7J is an approximation to the positive-definite discrete- 
velocity equilibrium ©. Same as with the isothermal 
models, this second-order approximation simply happens 
to be good enough for stable computations at Ma < 0.1. 
If we keep the relation p = (l/3)p to all the orders in 
j, we recover the exact positive-definite equilibrium of 
the isothermal 27-velocities entropic lattice Boltzmann 
model 0: 



fP( P ,j,p/3) = pWi 2 - ^1 + 3( Ja /p) 



a=l 
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Entropic stabilization procedure 0, 0, 0] can be ap- 
plied for the present model in order to achieve small val- 
ues of the transport coefficients but we do not address 
this here. 

Our approach to the discretization of the velocity space 
differs from the earlier considerations [2(il Elf . While 
we use the same Gauss-Hermite quadrature, we apply it 
on the grand canonical potential (J2J (that is, we eval- 
uate the velocity integral @ as pertinent to the mean- 
ing of a quadrature), and after that find the discrete- 
velocity equilibrium upon minimization of the discrete- 
velocit y gr and canonical potential (@J. Instead, authors 
of f2nL 12 ll| evaluate the local Maxwellian (that is, they 
evaluate a function, not an integral) of continuous kinetic 



theory, M{p,j,T;v) = p{2Tik B T/m) 



-D/2 



exp(— m(v 



u) 2 /2kBT), u = j/p, on the nodes of the quadrature. 
Certainly, just replacing M(p,j,T;v) -> Mi(p, j ,T;Vj) 
makes no sense because the conservation laws will be 
lost. However, it was noticed in |20j that when the 
second-order expansion of the Maxwellian is used instead, 
M^ 2 ' = pW(l + av a j a + bvaVfjjajfj), the replacement, 

+ bViaVifjjajp), coincides with 
the previously known second-order equilibrium of the 
isothermal lattice Boltzmann method [la, |l£| ■ It should 
be stressed that while in the Maxwellian must 

be truncated to second order (in order to rescue conser- 
vation laws), and thus the positivity of populations has 
to be sacrificed together with the second law of thermo- 
dynamics (Boltzmann's iJ-theorem), our equation (0 is 
just a good approximation to the positive equilibrium 10, 
and, if required, further terms can be computed in order 
to maintain positivity of the equilibrium populations to 
any required degree of accuracy. The discretization of 
the velocity space done at the level of generating func- 
tional (0) obviously violates none of the properties of the 
continuous kinetic theory. 

Ref. \2M indicated that the nodes of the fourth-order 



quadrature can be used for establishing a thermal model. 
Such model was indeed constructed and implemented in 
[l3l |2*^. However, even though the admissible temper- 
ature variation in this model is larger, it is bound to 
be less efficient than the lattice Boltzmann method. In 
that sense, the lattice Boltzmann model with energy con- 
servation derived in this paper is a good compromise for 
simulations of almost-isothermal low Mach number flows. 



In conclusion, we have given a microscopic derivation 
of a genuine lattice Boltzmann model with energy con- 
servation for simulations of incompressible and almost- 
incompressible flows. This model is a natural extension 
of the previously known isothermal models on the same 
lattices, and, due to the reasons explained above, it was 
entirely "overlooked" in all previous derivations. While 
we have considered here the 27-velocities lattice only, we 
shall address other important cases such as the 15- and 
19-velocities 3D lattices and the 9-velocity 2D lattice in 
our subsequent publications. The new models are as effi- 
cient as the previous isothermal lattice Boltzmann mod- 
els on the same lattices, and at the same time they extend 
considerably the domain of validity of lattice Boltzmann 
computations especially into the micro-flow domain. As 
we have already said, in principle, the spurious bulk vis- 
cosity of ILBM is not an obstacle for using them for in- 
compressible hydrodynamics simulation. However, even 
in that case, the present models can be preferred on the 
ground that they correspond more to the physics. The 
approach to the discretization of the microscopic theories 
uses explicitly the Legendre structure of thermodynamics 
[23j |. and can be used in other problems of reducing de- 
scription. This work was supported by the BFE-Project 
Nr. 100862. 
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